The Analytical Theory of Bulk Melting II: 
Variational Method Solution in the FCC Crystal 



O 

o 

Oh- 



Yajun Zhou * and Xiaofeng Jin t 
Surface Physics Laboratory & Department of Physics, Fudan University, Shanghai 200433, China 

(Dated: February 2, 2008) 

Continuing the arguments in Paper I (arXiv: cond-mat/0405487 ), we model the temperature de- 
pendence of interstitial defects in a surface-free face-centered-cubic (fee) elemental crystal and obtain 
the free energy and correlation behavior based on variational methods. We show that the avalanche 
of interstitial defects is the instability mechanism at the melting point that bridges Lindemann and 
Born criteria. 
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I. INTRODUCTION 

Bulk melting is a solid-liquid phase transition (SLPT) 
taking place in a surface-free crystal pj . The study of this 
phenomenon helps to clarify the essential driving force 
of the inhomogencous phase transition out of a homo- 
geneous system in that the surface-free condition rules 
out the influence of inhomogeneity at the surface bound- 
ary. Many endeavors have been made to test the previ- 
ous melting theories by probing into the behavior of bulk 
melting jj, especially in search of the relation between 
the widely-cited Lindemann 3] and Born |J| criteria. Lin- 
demann criterion proposes that melting is triggered by 
the avalanche of the root-mean-square (rms) atom dis- 
placement after it exceeds a threshold fraction (<5£) of 
the atom spacing (a) , where 5* L is called the critical Lin- 
demann ratio, a semi-empirical parameter once conceived 
as a lattice type characteristic; Born criterion argues that 
the vanishing of the shear modulus is responsible for the 
inability to resist lattice destruction at the melting point. 

In the recent molecular dynamics simulation con- 
tributed by Jin et al. |5|, it is demonstrated that the 
melting point of a surface-free Ar crystal does satisfy 
both Lindemann and Born criteria. In Ref. @, the Ar 
crystal melts when the shear moduli see a sudden down- 
fall (albeit not vanishing, which is consistent with the ex- 
perimental observation of residual shear modulus at the 
melting point @) and the atom displacement surges to- 
wards infinity simultaneously. Jin et al. have attributed 
this coincidence to the "Lindemann particle" (atom with 
displacement larger than 8* L times the atom spacing) clus- 
ters that emerge sporadically at low temperatures but 
abound throughout the crystal when melting point is ap- 
proached. It is found that within such clusters ("liq- 
uid nuclei" ) consisting of energetic "Lindemann parti- 
cles", the shear moduli difference is nearly vanishing: 
AOs = C44 - (Cn - C 12 )/2 re 0. Judging the role of 
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"Lindemann particle clusters" that satisfy both the Born 
and Lindemann instability criteria, it is then concluded 
in Ref. [5| that melting is governed by the "strongly- 
correlated" Lindemann and Born perspectives. However, 
several important questions remain unanswered by nu- 
merical results: 

Question 1: What kind of instability mechanism is di- 
rectly responsible for the non-zero shear moduli at the 
melting point? Is that contradictory to Born's original 
argument that melting is triggered by zero shear moduli? 
Question 2: How does a heterogeneous nucleation pro- 
cess manage to come out from a homogeneous system? 
Does the finite-size of the liquid nucleus at the melting 
point violate the principle "phase transition in the ther- 
modynamic limit" , which requires a system consisting of 
infinitely many atoms? 

Question 3: What are the most important factors, in 
terms of the parameters of the interatomic forces, that 
affect the magnitude of the parameter 8* L 1 How general 
is this S* L 0? 

In this paper serving as an extension of the model in 
Paper I to the three-dimensional (3D) face-centered-cubic 
(fee) elemental crystal, we report the detailed analyti- 
cal procedure that joins the Lindemann and Born cri- 
teria together in the SLPT and the mathematical argu- 
ments that address the three numbered questions above. 
We show that the spatial correlation between intersti- 
tial defects is the impetus that triggers and propagates 
instability in a crystal and that eventually undermines 
long-range order in the surface-free solid. Following the 
idea of defect-motivated transition in previous literature 
IEE3 [IHIHl and the previous understan ding s of defects' 
cooperation and aggregation [jj E, E, ll4L ITEl Hit ITU 

Eta, 

we use the J1-J2 lattice model plus vibrational 
Hamiltonian |2C| to formulate the cooperative creation 
and the spatial correlation of interstitials. Our results, 
which are based on fewer empirical parameters than that 
in Ref. not only epitomize the origin of instabil- 

ity and atom-scale pathway in the melting process with 
mathematical rigor but also help to elucidate the origin 
of "Lindemann particle clusters" and to explain the ex- 
perimental fact p| that 5* L is modulated by both lattice 
geometry and the profile of interactions. 

Paper II is organized as follows: Sect. II is a brief ac- 
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count of the methodology of our model for interstitial 
defects with a discussion of the symmetry of the Hamil- 
tonian; Sect. Ill uses variational method to work out a 
mean field approximation (MFA) solution to the three- 
dimensional (3D) model, where the temperature depen- 
dence of the concentration and the correlation of intersti- 
tial defects are investigated; Sect. IV discusses the phys- 
ical implications of the model solutions where analyti- 
cal results are compared to data from experiments and 
simulation; Sect. V summarizes the results in Paper II. 
Appendix A provides an alternative derivation of the for- 
mula which describes the temperature dependence of con- 
centration of defect; Appendix B gives the mathematical 
details involved in the evaluation of the correlation of 
defects. 



II. MODEL HAMILTONIAN AND ITS 
QUALITATIVE PROPERTIES 



A. Total Hamiltonian 



We begin our argument with the Hamiltonian describ- 
ing a system of N = 4£ 3 atoms, labeled as a — 1, . . . , N: 
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This Hamiltonian incorporates the kinetic energy of each 
atom labeled a (T a ) and pairwise potential energy be- 
tween atoms a and /3 (V a p), where V aa = 0, for a — 
1, . . . , N. For each atom labeled a, the potential energy 
which it experiences is the summation of the contribu- 
tion from all the other atoms, which reads Y^8=i Vaf3- 
From this sum, we extract a "vibrational potential" 
S^LiV^/?' which is harmonically oscillating with re- 
spect to interatomic distance when the atom a is per- 



turbed. The "configurational potential" 



formally defined as Ylp=i(VnP ~ Va/3 
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' vib ^ This Hamilto- 
nian separation method has been employed to study the 
lattice vibration's influence on the concentration of crys- 
tal defects 0. 

The model in Paper II basically treats an elemental 
crystal with fee structure (Ar, for instance). In order to 
take the interstitial defects into account, we use a lat- 
tice model with the standard NaCl-type structure. At 
absolute zero, the Na-like lattice is totally occupied by 
Ar atoms, while the Cl-like lattice is vacant and forms 



the totality of octahedral holes in a perfect fee crystal. 
Mathematically speaking, if one denotes a lattice posi- 
tion r by a triad (hkl), the Na-like (or Cl-like) lattice is 
characterized by odd (or even) h + k + I indices. At any 
finite temperatures, by disregarding the vibrations and 
lattice distortions (that is, by neglecting the atoms' ki- 
netic energy and its consequences) momentarily |22j , one 
may caricature the motion of the atoms in the solid as the 
stochastic hops on the two interpenetrating fee lattices. 
The configurations of the system are thus exhausted by 
all the possible ways to arrange N atoms of a fee elemen- 
tal crystal at 2N possible sites, so the summation over 
atom-pairs could be converted to a summation over site- 
pairs, as long as the site occupancy rate n at position 
r = (hkl) is included in the following way: 

a=l f3=l r r' 

1 U ( 

~ 2 E n hkl[Jl E 

h,k,l = l \ \h-h'\ + \k-k'\ + \l-l'\ = l 



nh'k'i' 



+J 2 



E 



n h n k , 



(2) 



/i"| + |fc-fc"| + |2-J"| 



Here, ^ r denotes summation over all the 2N sites, and 
n r denotes the number of the atoms occupying the site 
at position r, and n r — or 1. The "coupling constant" 
J™, nf in Eqn. © equals to V^ nf when the distance be- 
tween two atoms a and (3 is |r — r'|. We have used the 
following cutoff in Eqn. J™, nf = J\ when |r — r'| = 
the nearest-neighbor (NN) distance, namely, the distance 
between a pair of nearest "Na" and "CI" ; J r c °, nf = J 2 < 
when |r — r'| = the next-nearest-neighbor (NNN) dis- 
tance, namely, the distance between a pair of nearest 
"Na" sites (or "CI" sites); J™, nf = otherwise. In or- 
der that the NNN distance becomes the bond length in 
the fee crystal in our model (such as Ar, as opposed to 
NaCl), we require that J\ > J 2 and J 2 < 0. This J1-J2 
lattice model approximation is justified for interatomic 
forces that are both pairwise and short-ranged, which 
is physically applicable to solids where the interaction 
is governed by the Lennard- Jones (as in noble gases) or 
Morse functions (as in some metals) [T^. 

Therefore, in the fee crystal with 2N = 8£ 3 sites and 
N = 4£ 3 atoms (FIG.0, the configurational Hamiltonian 
now reads 
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rf = Eh* e 

h,k,l=l \ |/i-ft'| + |fc-fc'| + |i-H'l=l 



nh'k'i' 



+J2 n h „ k ,, vl \ (3) 

\h-h"\ + \k-k"\ + \l-l"\=2 / 

with the cyclic boundary condition: 

nh+2t,k,i — nh,k+2t,i = nh,k,i+M — n h,k,i (4) 
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FIG. 1: The way we label the lattice sites ( "*" ) and interstitial 
sites ("o"). 



that eliminates the surface and the atom number conser- 
vation condition: 



2/ 



(5) 



h,k,l=l 



that reduces the degree of independence for occupancy 
rate by one. Similar to the nomenclature in Paper I, we 
call the sites bearing the label (hkl) where h + k + I is 
odd (or even), namely, "Na" (or "CI") sites, as lattice 
(or interstitial) sites, or vice versa. The lattice sites 
and interstitial sites interpenetrate, as they did in the 
one-dimensional (ID) model in Paper I. 

In the 3D crystal, Ti vlh is taken into consideration in 
the form of vibrational free energy 



F 



vib 



3Nk B T log 



h(v) 



krfT ' 



(6) 



where k B is the Boltzmann constant, T is the absolute 
temperature, h is the Planck constant, (y) is the geomet- 
rical mean frequency of the crystal [2(1 |2lJ . It should be 
noticed that the formula above is valid in the Dulong- 
Petit limit - the temperature region where most melting 
processes take place and the heat capacity contributed 
by lattice vibration is asymptotically 3Nk B - 



B. Symmetry of the Configurational Hamiltonian 

We now study the symmetry properties of the config- 
urational Hamiltonian under the transformation: Ji i— > 
—J\- This will shed light on the properties of the "liquid 
phase" derived from our model Hamiltonian as well as 
the interplay of J\ and J 2 ■ 

First, by using the transformation uhh = ^ n hki — 1 
[23l IziL I23, one could map the model Hamiltonian in 
Eqn. (0 to the J1-J2 model for magnetism J2||, where 
o~hki = ±1: 



Q AT 

H conf (J 1 ,J 2 ,W)-^-(J 1 +2J 2 ) 



1 U ( 

h,k,l=l \ 



! ^ Oh'k'V 

h-h'\ + \k-k'\ + \l-l'\ = l 



+ J2 ^ ^h"k"l" I , 

\h-h"\ + \k-k"\ + \l-l"\=2 ) 

11 

0~h+2t,k,l = 0~h,k+2e,l = 0-h,k,l+2t = &h,k,U 

h,k.l=l 



0~h,kA = 0. 



(7) 

The partition function corresponding to this con- 
figurational Hamiltonian reads: Q conf ( Ji, J 2 , T) — 
£ w exp (-n coal (Ji, J 2 , {a}) /k B T) . 

Second, we notice that the Hamiltonian is invariant un- 
der any transformations that swap interstitial sites and 
lattice sites, i.e. mappings such as (hkl) 1— > (h+lkl), 
(hkl) i-> (hk + 11), (hkl) i-> (hkl + 1). (This is called 
"sublattice symmetry" 27].) It is conventional to as- 
sume |27j that there exists a temperature T' > OK, 
above which the partition functions Q conf (J%, J 2 , T) and 
Q conf (-Ji, J 2 ,T) both preserve the "sublattice symme- 
try" of the Hamiltonian. In other words, when T > T", 
both the (Ji, J 2 ) and (— J\, J 2 ) systems are characterized 
by equal amount of atoms occupying the lattice sites and 
interstitial sites. Since the number of atoms occupying ei- 
ther type of sites is N/2, we may infer that both partition 
functions (i.e. Q con{ (Ji, J 2 ,T) and Q conf (- J u J 2 , T)) 
are dominated by terms with the configuration {a} sat- 
isfying: 
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E 

h,k,l = l 



0~hkl 
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E 

h',k',l' = l 



oh'W = 0. (8) 



(-l)"-t-*-t-> =1 /_l\h' + k'+l' = _ 

So the transformation 



* I 1 \h+k+l 

0~hkl !-> °~hkl = V~ L ) a hkl 

sends one dominant configuration 



(9) 



M = y hk i,h,k,l£{l,2,...,2£} 

to another dominant configuration 
{a*}=\a* hkl ,h,k,l€{l,2,...,2£} 



21 



o-hki = 



h,k,l=l 



21 



(10) 



E < 



hkl 



h,k,l=l 



(11) 

in the partition function, and transforms the configura- 
tional Hamiltonian in the following way: 

H conf (Ji, J 2 ,a) ^ H cont (J 1 ,J2,a*) 

= H conf (- Ji, J 2 ,a) + WJi, (12) 
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where 



n conf ( Ji, j 2 , K}) - — ( Ji + 2j 2 ) 
1 2£ / 

g E °hfc2 ( ^1 E °h'k'i!' 

h,kj=l V Ift— /j'| + |fc-fc'|+U-J'|=l 



+ J 2 



E 



'h"k"l'' 



|/i-h"|+|fc-fc"|+|J-i"|=2 



(13) 



From the obvious identity: 

^-^ M conf (Jl,J 2 ,{g» 



E e 



■H conf (j 1 ,J 2 ,{ g *}) 

k BT ,(14) 



we see that when T > T", the configurational free energy 
(F conf = — /csTlog(2 conf ) has the following correspon- 
dence due to "sublattice symmetry" : 



gas model" |25j, that is, the properties of the system is 
predominantly defined by the NNN attractive bonding 
energy, and it is unnecessary to take the repulsion at NN 
distance into account. 

Therefore, we have obtained the "liquid phase" solu- 
tion to our model, which highlights the importance of J 2 
at high temperature. In the following section, we will 
show that for low temperatures T « f , our model be- 
haves differently as compared to the conclusions above, 
where different J\ parameter could change the properties 
of F conf (Ji, J 2 ,T) - 3N,h/2 drastically, indicating the 
loss of sublattice symmetry of the free energy. We will 
show that at low temperatures, the energy relationship 
J\ > J 2 < binds atoms together in the solid by break- 
ing the "sublattice symmetry"; near the melting point, 
the same energy relationship helps to propagate instabil- 
ity in the system and undermines the long-range order 
by revoking the "sublattice symmetry" . 
III. THE VARIATIONAL APPROACH TO THE 
3D FCC CRYSTAL 



37V 

F co » i (J 1 ,J 2 ,T)- — (J 1 + 2J 2 ) 
= F conf (- J u J 2 , T) — -— {-J x + 2J 2 ) . (15) 



This equation is parallel to Eqn. (12) in Paper I, which is 
a formula invariant under the transformation J\ i— » —J\. 

The analysis above reveals that when T > T', the free 
energy of the system (up to a ground energy term that 
reads (3iV/2)(±Ji + 2J 2 ) ) is sensitive to J 2 but proba- 
bly insensitive to J\ - it is because F coni (±J 1 , J 2 , T) — 
(3iV/2) (± Ji + 2 J 2 ) remains unchanged even when Ji al- 
ters its sign, manifesting the irrelevance of J\ in this 
temperature region. This scenario is consistent with the 
previous description of the liquid phase in the "lattice 



A. Phenomenological Parameter Expansion of the 
Free Energy Functional and the Model of Melting 

In our probe into the 3D fee elemental crystal, we aim 
to obtain a solution based on modified MFA. We first 
define a random variable called the local order parameter 



L(v) = (-l) h+k+l a hkl 



( _ 1)h+k+l {2nhki 



1) 



(16) 



where site r bears the integer label (hkl). According to 
the Landau-Ginzburg theory of MFA [2^] , the free energy 
functional contributed by configurational Hamiltonian of 
the 3D fee model 7i conf could be expanded phenomeno- 
logically as 



^ COIlf [L (r)] = -TS coni [L (r)] + \ (^j J d 3 r { [(l - I? (r)) Hi + (l- L 2 (r)) 2 H 2 ] + 7 [VL (r)] 2 } 

(17) 

where (V2/a) 3 J d 3 r =X) r denotes the summation over all the 2N sites. (Recall that the NN distance in the 2N 
sites is a/y/2, so the volume of the smallest cube on the lattice is (a/v^) 3 , which forms the ratio between the site 
summation J2 r an d tne spatial integration J d 3 r.) Hi (> 0) and H2 (< 0) are two energy parameters to be determined 
later as functions of J\ and J 2 . 7 (> 0) is the phenomenological "domain wall energy" coefficient. The configurational 
entropy S coni [L (r)], which appears in the equation above, reads as the entropy of mixing pll |24|: 

3 



S cont [L (r)] = -fc fl 




(18) 



F \L (r)] is contributed by vibrations with geometric mean frequencies (f;) (lattice mode) and (z/,-) (interstitial 



mode) 



which reads 



[L(r)} = 



3k B T ( \/2 



d 3 r 



L 2 (r)) 



log- 



const 



(19) 
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when expanded in power series of (l — L 2 (r)J. 

Now the total phenomenological free energy functional F [L (r)] reads: 



F [L (r)] = F conf [L (r)] + F vih [L (r)] = (^j J d 3 r / (L 



(r),VL(r)) 



(20) 



where 



l r 



/(L(r),Vi(r)) = - (l-L 2 (r))H 1 + (l-L 2 (r))^ 2 -3A :B T(l-L 2 (r))T 



2 



l - L ^ log l^) + l±M log l±^ 



| |VI(r)] 2 



(21) 



and T = (1/4) log {(vi) / (vi)) > is a factor modeling 
the "lattice softening" effect, that is, interstitials have 
a lower vibration frequency than normal atoms at the 
lattice sites. By variational methods, we may apply the 
Euler-Lagrange equation 



3/(L(r),VL(r)) „ Of (L (r) , VL (r)) 



dL (r) 



dVL (r) 



(22) 



to optimize the free energy functional with the boundary- 
free condition and obtain the equation of motion for L (r) 
in the "Poisson equation" form: 



p(L(r)) = - 7 V 2 L(r), 



(23) 



where 



p(L(r)) = HiL (r) - 2H 2 L (r) (L 2 (r) - l) 

-3k B TL (r) T - fc B Ttanh _1 L (r) . (24) 

The "phenomenological electric charge density" 
p(L(r)) immediately gives rise to a simple mathemat- 
ical model of the catastrophe of the long-range order: 
For sufficiently low temperature T, the curve p(L(r)) 
intersects the positive L (r)-axis at least once (FIG. [2J , 
so there is a homogeneous distribution L (r) = (L)t ^ 
that makes p(L(r)) vanish. When T — T m , where T m 
suffices 



p((L) q 



= and 



dp(L(r)) 



dL (r) 



= 



L(r) = (L) 1 



simultaneously, the intersection becomes a tangent point 
and we have the following equations: 



3T 



(L) T 



~Hi-2H 2 (3{L) 2 Tm -l) 



k B T„. 



- tanh 1 (L)r 



4H 2 
k B T„, 



(25) 



In FIG. El k B T m = Ki/2.53. For temperatures higher 
than T m , p (L (r)) no longer intersects the positive L (r)- 
axis, suggesting the loss of long-range order. From 




FIG. 2: (color online) Local properties of p(L(r))/fcsT 
near its node: {L) T . In this specific case, Hi = 6 1 7-^2 1 , 
T = log (4/3) . The tangent node {{L) T , 0) is related to the 
melting point. 



Eqn. l(2*oT) we know that as H. 2 tends to zero, so does 
(L) T . Therefore, in the absence of "virtual attraction 
between defects" (oc H.2, to be elaborated immediately 
in next two subsections), first-order SLPT is not possible 
according to this model. 



B. The Cooperation Effect and the Ti.2 Term 

In the MFA scenario, 
F[L(r)} 

1 / Fx 

7>vib 71 nconf ill i \ A$ 



+ / d 6 r {(6J1 (n r rv)) 



+UJ 2 (n r n r „) + 7 [VL(r)] 2 }, 



(26) 



where n r = (<r r + 1) /2, r and r' are NN sites, r and 
r" are NNN sites. One simple-minded approximation 
yields the result: (a r a r i) = — (a r a r ") = —L 2 (r). This 
approximation is equivalent to the following statement: 
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when the correlation between sites is negligible and L (r) 
distribution is uniform, the lattice (or interstitial) site 
occupancy is (1 + \L (r)|) /2 (or (1 - \L (r)|) /2), so that 
multiplication rule in probability theory infers that |24| . 



(n r n T i) 

l + |L(r)|l-|L(r)| 1 — L 2 (r) 



(27) 



2 (( n r"r")lattico sites + ("r«r" ) interstitial sites) 



1 + \L(t) 



1-\L(t) 



1 _ 1 - L 2 (r) 

2 4 ' 



(28) 



Such an argument is reasonable when \L (r)| is sufficiently 
close to 1 and the interstitial concentration is low enough 
to obscure the cooperation between defects. When T = 
0, the approximation above gives the defect concentration 
in the form of (1 - |£(r)|) /2 ~ exp (-u/2k B T) where 
u = 6 Ji — 12 J 2 is the excitation energy of one interstitial 



defect. This asymptotic behavior of the partition func- 
tion agrees with the well-established theory of Frenkel 
defects where excitations of defects are assumed to be 
spatially independent pcj . 



However, the approximation (cr T a r 



(a r a r ") is far 



from accurate when the interstitial concentration is high 
and correlation between atoms should be taken with care. 
By the qualitative analysis of the exact partition func- 
tion, we see from Eqn. JSJ) that (cr r <T r /) vanishes in MFA 
scenario when T > T', so that F conf (J 1 ,J 2 ,T) - ZNJx/2 
is independent from J\ when T > T' while employ- 
ing MFA. However, it does not follow that (a T a r ") = 
— (c r ov) = 0. We thus have to make corrections to the 
estimate of (cr r o" r "), especially when |£(r)| is far away 
from 1. To see the indispensability of this correction, 
we re-examine the partition function from the graph- 
theoretic perspective. In the thermodynamic limit, when 
T > T", by replacing the summation over all states with 
the summation over thermodynamically most probable 
states (which is a "steepest-descent argument" similar to 
Eqn. (6) in Paper I), we may write the exact partition 
function as: 



lf (Ji, J 2 ,T) = ^exp 
M 

(Ji + 2J 2 ) 



H c 



(Ji, J 2 ,ct) 



k B T 



X) exp 

{a hk l=±l} 



(Jl, J2,er) 



k B T 



exp 



2k B T 



Eq r * tanh r — tanh 



cosh 2Ar -^3- cosh 2Ar 

LObn 4k B T LOhn 4k B T r,s 



2 2N 



exp 
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" 2k B T 



(Ji + 2 J 2 ) 



Lusn 4k B T 



^ 9o,a tanh J 



4fc s T 

J 2 
4fc B T 



4fc B T 



(29) 



where g r ^ s is the number of loops that contain r antibonds (lines that join two NN sites) and s bonds (lines that join 
two NNN sites) and r must be an even number because every loop is closed. (The evenness of the number r confirms 
again the Ji i— ► — Ji symmetry, cf. Appendix A in Paper I) The final "~" in the equation above results from the fact 
that tanh( Ji/4fcsT) i=s 0, and that cosh( Ji/4ksT) i=s 1 since the temperature T is high enough. It follows that 



J ,cont (Ji,J 2 ,T) 



3A 



( Ji + 2J 2 ) 
3A 



-fc B TlogQ cont ( Ji, J 2 ,T) - — ( J x + 2J 2 ) 



= 2Nk B T log cosh -p— - 2Nk B T\og 2 - k B T log ( V g , s tanh s 
4k B T 



4k B T 



(30) 



still varies with respect to T when T > T', and this 
variation is dependent on the NNN attractive interaction 
J 2 , reinforcing that (<J r o~ r >>) is still non-zero when T > T' . 

From the behavior of the exact partition function, we 
see that the mean NNN correlation is more persistent 
than its NN counterpart. In other words, the system 



tends to "preserve" more bonds than assumed in the ap- 
proximation (<7 r <7 r <) = — (<o~ T <J T »). The inaccurate ap- 
proximation (o~ r o~ r >) = — (crov) infers that when intcr- 
stitials are excited, the creation of 6 antibonds is accom- 
panied by annihilation of exactly 12 bonds, which is not 
necessarily the case. 
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In order to offset the overestimate in the bond annihi- 
lation, we have to write down 

6J1 (n T n r >) + I2J2 (n r n T ») 

= 6J 2 + (l-L 2 {r))H 1 + (l-L 2 {r)) 2 H 2 (31) 

where Hi = (3/2) (— 2J 2 + Ji), and the negative coeffi- 
cient Ti 2 , which is proportional to J 2 , is to be determined 
in the next subsection. 



The Phenomenological Parameters H2 and 7 as 
Functions of Ji and J2 



sites are scarcely occupied by two atoms simultaneously, 
NN correlations are immaterial as compared to NNN cor- 
relations. In the light of this, we may retain Eqn. 1)27(1 in 
the final expression of free energy functional. Comparing 
this with Eqn. I|3I|I . we have 



n 2 = -J2 < 0. 



(34) 



When spatial inhomogeneity is significant, a domain- 
wall energy proportional to [VL (r)] 2 should be taken into 
account in order to offset the miscalculation (A) based 
on the "mean- field" assumption L (r) = L (r') = L(r"). 
In Eqn. 1(21)1) . a good estimate of 7 is a 2 Hi/2 because 



When the concentration of interstitials is considerable 
(about 3% ~ 5%), the creation of di-interstitials (exci- 
tation of two interstitial atoms at NNN distance) helps 
to preserve more NNN atom pairs than predicted in 
Eqn. (|28|l . This cooperation effect inherent in the Hamil- 
tonian lowers the energy cost in exciting interstitials be- 
cause a di-interstitial has less formation energy than that 
of two interstitials at farther separate distance in our 
model 30 . Accordingly, we have to correct Eqn. I|28|l 
by adding a positive term proportional to (1 — L 2 (r)) 2 
(symmetry and analyticity of the free energy preclude 
terms including L (r) or |L (r)|). One scheme to outline 
this correction is by estimating the percentage number of 
occurrence of NNN-contact holes in the lattice sites as 



12 
~2 



1 



\L{v) 



\^-L 2 {r)) 2 



(32) 



where "12" is the coordination number of NNN pairs. 
FIG. 01 shows how "virtual attraction" [£j between de- 
fects is taken into consideration based on probabilistic 
arguments. The system preserves more bonds than in 
the simple-minded approximation: when two holes on 
the lattice sites are close (in NNN contact) instead of 
being farther separate, they save the bond energy in the 
amount of J 2 . In total, the saved bonds by stochastic 
contact of holes lower the system energy in magnitude of 
(3/8)|J 2 Kl-L 2 (r)) 2 . 

Dynamically speaking, this correction could be also 
understood by the NNN correlation on the lattice sites. 
Since the occupancy of lattice sites is approximately 1, 
the relaxation time of bonds on lattice sites is much 
longer than that of the bonds on the interstitial sites 
and the antibonds. So the dissociation of bonds is 
a "slower reaction" than the creation of antibonds. 
That is why (^r^V') lattice sites should be greater than 
[(1 + |i(r)|)/2] 2 . " ' " 

Therefore, after considering defect attraction or NNN 
correlation on the lattice sites, the ensemble average 
(n r n r ") should be modified into 



(n r n r //) 



i^ + §(l-l'M)'. (33) 



so as to account for the crystal's tendency to retain as 
many NNN atom pairs (bonds) as possible. Since NN 



A(2Ji (<j t <t t .)) 



2 I p 

- 1 L(r) L(r') 



i(r) i(r>; 



£(r) L(r') 



- J\ (2pr — p 2 — r 2 ) 



J 1 {p-rY = J 1 -[VL{v)Y 



A(2J 2 {(J r (T r n)) 

£(r") \ / i(r") 

J, y \ I J, 



(35) 






- J 2 a 2 [VL (r) 



(36) 



The total domain-wall energy contributed by the non- 
vanishing VL (r) at site r is thus 



1 a 1 t a 
- x 6 x - Ji — 

4 2 V 2 

yWi [VL(r)] 2 . 



J 2 a 2 



[VL(r) 



(37) 



Here, the factor 1/4 results from the transformation 
cr r = 2n r — 1, the factor 1/2 is applied to offset repeated 
counting, the number 6 is a consequence of coordina- 
tion number and the Pythagorean theorem applied to 
the VL (r) vector decomposition. 

Physically speaking, at temperatures far lower than 
the Dulong-Petit limit, the domain wall energy could be 
interpreted as the energy contributed by the dislocations 
which destroy the continuity of L (r) on the interface. 
The (7/2)[VL(r)] term associated with rare disloca- 
tions contribute little to the total free energy at low tem- 
peratures. When temperature is sufficiently high, large 
fluctuations of L (r) will be commonplace and the gradi- 
ent term will become indispensable if we want to evaluate 
the total free energy correctly. (It should be emphasized 
that the vibration term containing T arises from entropy 
effect and does not contribute to domain wall energy. 
Even when the distribution of L (r) is inhomogeneous, 
the average vibrational energy for every degree of free- 
dom is the same, say, ksT, regardless of vibrational fre- 
quency. Therefore, the expression of 7 is not dependent 
on T.) 
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(a) 



in a dilute "gas" of interstitial monomers and intersti- 
tial oligomers, in search of a minimal free energy cor- 
responding to an optimized long-range order parameter 
(L) T (FIG.S(a)). 

Second, we notice that the aforementioned "colli- 
sion process" causes a fluctuation of L (r) , governed by 
Eqn. (|23|) that resembles the equation for a globally neu- 
tral plasma when T <C T m . (This is because p(L(r)) 
changes sign as L (r) varies in the vicinity of (L)t, which 
is analogous to the coexistence of positive and negative 
charges in a plasma.) The Green function of fluctuation 
response G(r, r') satisfies 




(b) 




V2 



FIG. 3: (color online) (a) This shows two layers in a crystal. 
In the "upper layer", two holes are separate whereas in the 
"lower layer", two holes are close, (b) This shows the cross- 
section view of the "upper layer", with loss of eight bonds 
(in dashed lines) as compared to a perfect layer, (c) This 
shows cross-section view of the "lower layer", with loss of 
seven bonds (in dashed lines) as compared to a perfect layer. 
We could see that two close holes save a bonding energy in 
the amount of J2 (bonds are denoted by thick black lines in 
(b) and (c)). 



D. "Chemical Equilibrium" and its Demise Due to 
Locally Isotropic Instability 

The physical implications behind the mathematical 
model outlined in Sect. III. A provide much more infor- 
mation about the melting pathway. 

First, for T < T m , we find that p ((L) T ) = is equiva- 
lent to the "chemical equilibrium" condition: 



Here. 



1 — [defective cell] 

1 + |(L) T | [non-defective cell] 

' denotes equilibrium concentration, and 



(38) 



As = 2 {Hi - 3k B TT) \(L) T \ - 4W 2 |(L) T | {\(L) T \ 2 - 1) 

(39) 

denotes the energy difference between two types of crys- 
tal cells: the defective cell incorporating (in the statisti- 
cal parlance) more than one (inclusive) interstitial and 
the non- defective cell including less than one intersti- 
tial (inset of FIG. 0] (a)). (In Appendix A, we will de- 
rive the form of As from a perspective independent from 
the arguments in previous subsections.) This dynamical 
equilibrium is achieved by stochastic creation and anni- 
hilation of interstitials, or vividly speaking, the collisions 



7V 2 



dp(L(r)) 



dL (r) 



G(r,r') = -fc s T5(r-r'), 
(40) 

where S (r — r') is the Dirac delta function (See Ap- 
pendix B for the derivation of the Green function). 
For sufficiently low temperature T <C T m , the partial 
derivative in the equation above is always negative, so 
G (r, r') is a propagator that decays exponentially to 
guarantee the stability of the homogeneous distribution 
L (r) = (L) T . However, as T approaches T m from be- 
low, it is possible that dp(L(r)) /dL(r) > for certain 
values of L(r) not far from (L) T , which results from a 
non-Gaussian perturbation (Appendix B) to the homoge- 
neous distribution. In the complex wave-number plane, 
apart from the poles at purely imaginary wave-numbers 
(decay mode), the Fourier transform of the Green func- 
tion would then also encounter real wave-number singu- 
larities (oscillation mode), indicating that homogeneity is 
only preserved for a finite volume within which fluctua- 
tion propagates as a sinusoidal wave. When T = T m , 
which is defined in Eqn. (|25[) . the "phenomenological 
electric charge density" p(L(r)) no longer changes its 




.0 0.1 0.2 0.3 0.4 0.5 



0.75 0.80 0.85 0.90 0.95 

(a) 



4.5 0.6 0.7 0.6 0.9 1.0 
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



(b) 



FIG. 4: (color online) (a) This plots (L) T — T dependence 
for the same TH.x/\H%\ and T as in FIG.|2] where catastrophe 
occurs at the "melting point" T m — 0.89| J2 |/fc.B. Inset shows 
a "chemical equilibrium" between non-defective (up) and de- 
fective (down, interstitial defect highlighted by red) cells at 
T < T m . (b) Plotted here is the dimensionless parameter 
pi/kBT m as a function of {L) T (inset shows details of the 
function using different scales), pz appears in the expression 
of 7r£, the critical size of the liquid nucleus. 
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sign in the vicinity of (L) T , and could be expanded as 
p (L (r)) = — p? (r) pi — p? (r) p 3 + ■ ■ ■ in the vicinity of 
(L) T , where p (r) = L (r) — (L) T . 

Physically speaking, this expansion of p (L (r) ) infers 
that the uniform distribution L (r) = (L) T is unstable 
when there is a radial symmetric perturbation (loo wave) 
of L (r) judging the Gauss' theorem: 



S=dV 



VL (r) • dS 



p(L (r))dV 



> 0, L (r) > 
<0, L(r)<0 



(41) 



When a perturbative Y 00 wave initiates, the surface in- 
tegral boosts when the spherical volume is increased, so 
the value of L (r) grows monotonously along the radial 
direction. Anisotropic Yj m modes do not apparently lead 
to steady growth of L (r) values and the related propa- 
gation of instability. It is because for Yi m modes where 
I > 0, there is no such a causal relationship between the 
surface integral's growth and the growth of the L (r) in 
the radial direction. 

It can be verified mathematically that anisotropic Yi m 
modes are quenched as ~ r l in the short range and do not 
account for the instability. This is because p (r) satisfies: 



7VV(r) =P 2 (r) P2 + P 3 (r)p3. 



(42) 



For anisotropic excitations, where the asymptotic behav- 
ior is p(r) ~ p (r) Yi m (6, <fr), r — |r| — > (the original 
point of r is placed at the center of excitation), I > 0, we 
find 



7X "(r) 



/(f + l)7X(r) , P2X 2 (r)Y h , 



r 



+ 



^ 1(1 + l) 7X (r) 

where x ( r ) = | r | A* ( r )- 

p (r) ~ r l — > 0, as r — » 



(43) 



(44) 



Meanwhile, Yoo wave has the following asymptotic be- 
havior: 



p{r) — ^ 0, as r -> 

P3 



(45) 



Therefore, isotropic excitation dominates the instability 
mechanism at T m , because all the Y; m (I > 0) modes are 
overwhelmed by the Yqo mode in the short range. 



E. Relations with Born and Lindemann Criteria: 
"Lindemann Particle Clusters" Revisited 

In the context of instability induced by locally isotropic 
excitations, "quasi-neutrality" in the "plasma" is at- 
tained by establishing L (r) > and L (r) < domains, 



each in diameter of 7r£ ~ \[2ali\Tr'(2pj > kBT m )~ 1 (Ap- 
pendix B elaborates on this estimate, and FIG.QIb) plots 
the dimensionless parameter p^/ksTm), which is the av- 
erage size of the locally isotropic excitations of atom clus- 
ters. These highly cooperative and energetic atom clus- 
ters gives rise to a catastrophe of global long-range order 
at Tm'- 



\(Lh 



0. 



T > T m 







(46) 



The vanishing long-range order at temperatures 

higher than T m casts the system into sublattice symme- 
try. 

Such locally isotropic excitations near T m form spheri- 
cal domains of instability (SDIs) and amount to one pos- 
sible interpretation for the origin of "Lindemann parti- 
cles" |5j and their aggregation. In such SDIs, due to 
the nature of the fluctuation and the relaxation, the dis- 
placed atoms are moving collaboratively, energetically 
and isotropically, so that we could regard them as the 
physical entity that bridges Born and Lindemann criteria 
at the melting point according to the following argument: 
On one hand, the totality of atoms in each SDI exhibits 
isotropy, which guarantees exactly vanishing shear mod- 
uli difference: ACs = - that is, SDIs satisfy Born 
criterion; On the other hand, the displacement of the 
atoms in all the SDIs exceeds a8* L - that is, the average 
displacement of all atoms in the Euclidean space, which 
packed spheres cannot perfectly fill. In the light of this, 
the atoms in SDIs are "Lindemann particles" by definition 
- that is, SDIs also satisfy Lindemann criterion. 

In short, the locally isotropic plasma instability at 
the superheating limit ( "melting point" of a surface-free 
solid) T m is initiated and propagated by SDIs, and T m is 
exactly the temperature at which Born and Lindemann 
criteria coincide. These SDIs declare the demise of the 
chemical equilibrium originally established by collisions 
in a dilute "gas" of interstitial defects, and undermine 
the long-range order in a solid in the mean time. 



IV. DISCUSSIONS 

From the arguments in the last section, we have 
provided accounts for the instability mechanism (why), 
molecular pathway (how) and phase transition temper- 
ature (when) related to the melting of a surface free fee 
elemental crystal. The answers to why and how also help 
to clarify the physical basis for the equivalence of Born 
and Lindemann criteria. This section will provide further 
implications of our theoretical model and give answers to 
the numbered questions in the Introduction. 

The crucial instability mechanism (why) of melting is 
found to be the cooperative creation of interstitials that 
caused the avalanche of displaced atoms. The coopera- 
tion effect transfers the information of one displaced atom 
isotropically to its neighborhood and finally results in the 
formation of SDIs. The residual shear modulus at the 
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melting point is due to the inability of spherical domains 
to fully occupy the Euclidean space. Born's argument 
is still valid within each SDI, but no longer valid in the 
melting crystal as a whole. This answers Question 1. 

Therefore, for the whole crystal, lattice softening or 
vanishing shear modulus is not as decisive in the melting 
process as the cooperation effect is. In the melting point 
formula Eqn. (|25|) , it is clear that the first order transition 
disappears when TI2 is zero. The cooperation between 
the interstitials proves to lower the energy cost to excite 
defects thereby paving the way for destroying the long 
range order in the lattice. For fixed bonding energy J2 
and lattice softening T, T m is lowered when \H.2\ / 



ratio is enhanced, which is in accordance with the ID 
exact solution which relates high | J2 1 / | J\ | ratio to great 
disorder. 

To justify this, we will prove the following 
Theorem: 



dT„ 



> 



1 / h 2 ,t 



(47) 



in the 3D fee case. 

Proof: From Eqn. (|25|) . pick (L) T > 0, we may arrive 
at the following conclusion: 



d 



Hi - 3k B TT 



d(L) Tm V \H 3 
Therefore, we have the following argument 



4 <£>T m ( 


1 - < L >L) 


2 

tanh (L) T 




i 1 wk] 


1 2 

tanh" 1 (L) Tm 



{L) 7 



8x 4 da; 
(1 - x 2 Y 



< 0. 



(48) 



\H 2 



is fixed 



Eqn 



^Ti x - 3k B TT / 
|T| is fixed, T m / 



Hi 



Thus, (dT m /dHi) H2>r > is evident. ■ 

That is to say, T m is not only determined by the 
bonding energy related to H.2, but is also modulated by 
the difficulty to create interstitials, which is signified by 
Ti.i - the energy barrier that hinders information trans- 
fer between atoms. The mathematical inequality above 
could be physically interpreted as "the better information 
transfer, the lower melting point" . Unlike the liquid-gas 
phase transition involving the thorough dissociation of 
bonding atom pairs, in which J2 plays a decisive role 
in determining the boiling point, the SLPT is possible 
when destruction of the lattice structure and the excita- 
tion of interstitials are both highly encouraged. There- 
fore, some materials could have high boiling points and 
"disproportionately" low melting points, if they have a 
large and a small \Hi\- (This may have shed light 
on the peculiar behavior of Ga and In, two metals that 
melt near room temperature and boil at thousands of 
Kelvins, although neither metal falls into the category of 
fee structure.) The lattice softening T > is also con- 
ducive to melting in that Ti\ — 3k B TY < TL\ enhances 
the apparent \H,2\ / \Hi\ ratio. From this we know that 
the atom displacement and lattice softening are mutu- 
ally complementary and this rule underlies the intrinsic 
relation between Lindemann and Born criteria. 

The atom-scale pathway (how) of melting is spa- 
tial correlations of atom occupancy fluctuations. The 
correlation-response scheme proves to be the effective 



way to transfer information in the 3D system discussed 
in this paper. The anomalous sinusoidal correlation wave 
results from the non-Gaussian fluctuations of atoms. In 
Ref. p| , it was observed that melting is preceded by the 
non-Gaussian behavior of atom displacement, but the 
causal relationship between the deviation from normal 
distribution and the SLPT was not explicitly established. 
This paper presents a clear demonstration of this causal- 
ity in Sect. III.D and E. From the conclusions in those 
subsections, we may find the critical diameter of the SDIs 
for the example shown in FIG. (b): 7r£ m i„ = 15.03a, 
and the spherical volume 7r(7r£ m i n ) 3 /6 = 1780a 3 which 
contains 890 atoms in average. This forms a good com- 
parison with the statement in Ref. [f| that the "Linde- 
mann particle" clusters consist of 10 2 — 10 3 atoms. The 
critical volume of the liquid nucleus 7r(7r£ m i n ) 3 /6 will be 
evaluated in detail in Appendix B. 

By exploring the atom-scale pathway in the thermody- 
namic limit, we highlight the fluctuations that propagate 
as a sinusoidal wave. This wave cuts the space into sep- 
arate compartments exhibiting ostentatiously mutual in- 
dependence, and is thus responsible for the illusion that 
the inhomogeneous instability suddenly emerges from 
nowhere. There is no contradiction between the nuclc- 
ation and the thermodynamic limit, after all. This an- 
swers Question 2. 

The melting point (when) is defined in Eqn. I|25|) and it 
is shown to be the superheating limit predicted by Born 
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and Lindemann criteria simultaneously. When rewritten, 
the melting formula takes the form 



3\J 2 \{L)\ 



< 



9|J 2 



The inequality results from the following 



(49) 



Theorem: 



(L) Tm 
2 1-(L)% 

3< 1 



tanh -1 (L)r 



(50) 



Proof: Pick (L) T > 0, 



1 - (L)j 



tanh- 1 (i) Tm - | (L) 3 Tm = J 



^T m 2x 4 (2 - x 2 ) dx 
(l-x 2 ) 2 



> 



(51) 



We can check the reasonability of the formula above by the following table: 





Ne 


Ar 


Kr Xe 


Rn 


Cu 


Ag 


Au 




/K 




24.56 


83.8 


115.79 161.4 


202 


1357.6 1234.93 1337.33 


T h 


/K 




27.07 


87.3 


119.93 165.1 


211.3 


2840 


2435 


3129 


A-ffvap 


/kJ mol 


-l 


1.7326 


6.447 


9.029 12.636 


16.4 


300.3 


250.58 


334.4 


AH lus 


/kJ mol 


-l 


0.3317 


1.188 


1.638 2.297 


2.89 


13.05 


11.3 


12.55 


ASf us 


/J mol - 


X K -1 


13.5 


14.2 


14.1 14.2 


14.3 


9.61 


9.15 


9.38 


\J2\N A 


/kJ mol 


-l 


0.3441 


1.273 


1.778 2.489 


3.215 


52.23 


43.65 


57.83 


9|J 2 | 
4fe B 


/K 




93.109 344.37 481.13 673.58 870.07 


14134 


11812 


15649 


T 


/K 




29.47 


100.6 


138.9 193.7 


242.4 


1629 


1482 


1604 


1(1-1 


< £ >tJ) 




0.0475 0.0440 0.0435 0.0430 0.0420 


0.0180 


0.0195 


0.0160 


AS m 


/J mol - 


i K -i 


8.35 


8.53 


8.56 8.57 


8.63 


10.0 


9.93 


10.2 



In this table, all the elements (Ne, Ar, Kr, Xe and Rn from the group; Cu, Ag and Au from the IB group) 
assume fee structures in the solid phase. The first four rows were taken from the descriptions of individual elements 
in the online encyclopedia, http://en.wikipedia.org/ (accessed Sept. 10, 2004). Tg is the equilibrium melting point 
of crystals with surfaces. Xb is the boiling point of the liquid. A_ff vap is the heat of evaporation. AiJf us is the heat of 
fusion. The data in the other rows are obtained/estimated as follows: ASf us = AH[ us /Te is the entropy change due 
to fusion; T m — 1.2Tg is an estimate of superheating limit based on Ref. [5j; | J2I = (AiJ V ap + A_fff us )/6A^4 is a rough 
estimate for the J2 parameter, where Na is the Avogadro's number; |(1 — \ {L) T |) is the critical concentration of 
interstitial defects that initiates melting, where \ {L) T | is calculated from | J2 1 , T m and Eqn. I|25l) : 



AS„ 



~»conf I 



\L(r)=0 



conf I 



= N A ki 



L(r)=\(L) T J 

21og2+(l- |<L) T J)log 



In principle, 



+ (1 + I (L) Tm I ) log — = 



(52) 



(53) 



where V and V s are the molar volume of liquid and 
solid, respectively. Therefore, the disagreement between 
ASf us and AS m may be explained as the volume change 
and the T m /TE ratio that are not fairly covered in the 



above calculation of AS m . In general, the table justifies 
the inequality T m < 9 | J2I /4fcs, and shows similar ASf us 
(or AS m ) values for the elements in the same group. 
The data also suggest that melting is triggered when the 
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concentration of interstitial defects is very low but still 
high enough |3l| to induce correlation and cooperation. 

For noble gases Ne, Ar, Kr, Xe, Rn, the interac- 
tion mode is the well-known van der Waals force and 
is described by the Lennard- Jones function |5j, so they 
share theoretically the same (and practically similar) di- 
mensionless parameters \Hi\ / \H\\, T and V 1 /V s in our 
model. Hence they should have the same dimensionless 
values: \{L) T |, \J2\/k B T m and AS m /NAk B . Another 
dimensionless parameter 5* L could be obtained by using 
the following 
Theorem: 



ik B T m = A^{v I )m{a8l) 



(54) 



where (v 2 ) is the average frequency 0. 
Proof: By multiplying the Newtonian equation mr(t) — 
F (i) with r(t) and taking time average, we can establish 
the following Langevin equation 



mr(t) ■ r(t) = r(t) •F(i), 



(55) 



where m is the mass of the atom, on which a force F (t) 
is exerted at time t. 



r(t) • r(t) 



lim — 

r—*oo T 

lim — 

r— >oc T 



r(i) • ¥(t)dt 



r(i) ■ df(i) 



1 1 r 

= lim — [r(t) • r(i)ir =n ] — lim - / i 2 {t)dt 

r— >oo 7 t^oo t Jq 

(56) 

r(i) • i"(i)|J" =0 is finite insofar as vibration and Brownian 
motion are concerned. So lim T ^ 00 (l/r) r(t) ■ r(i)|[ =0 = 
0, and r(t) -r(t) = —r 2 (t). According to the ergodic- 
ity argument, we may replace the time average by the 
ensemble average as long as T < T m . Therefore, 



mv(t)-r(t) = -mr 2 (t) = -m(r 2 ) = -3k B T 
II 

r(t)-F(t) = (r ■ F) = -4tt 2 to <V 2 r ■ r) (57) 
Define {v 2 ) = (^ 2 r 2 ) / (r 2 ) Q, then we have 

lim 3k B T = lim 4tt 2 to (v 2 ) (r 2 ) 



3k B T m = 4TT 2 m(v 2 ) (a5* L y 



(58) 



For Lennard-Jones potential, m(y 2 } oc |7^2 1 oc T m , so 
there is a universal 5£ for all fee crystals governed by van 
der Waals forces. In our model, k B T m j \ J 2 \ is a function 
of J (L) T J that varies slowly with respect to | (L) T | . The 
variation itself infers that 8* L is not a universal constant 
independent of interaction force, and the slow variation 
explains the reason why similar materials have nearly 



identical <5|/s. In a nutshell, the most influential factors 
for the critical Lindemann ratio 5* L are the profile {not 
every detail!) of the interatomic force that determines 
(L)T m as a function of Ji, J2 and T in different crystals. 
There is no such a strictly universal <5£ for fee lattice • 
This answers Question 3. 

We admit that the current model solution may need 
further refinement in terms of parameter estimate. The 
relations H2 = (3/8) J2, 7 = a 2 H.i/2 and the representa- 
tion of 7r£ min in this work are obtained by bold assump- 
tions and simplifications. However, judging the physical 
origin of the collaborative effect and domain wall energy, 
it is clear that H2 < 0, 7 > is not questionable even 
when they are evaluated exactly. p 3 > is also guaran- 
teed in all cases: 
Theorem: p% > at T = T m 
Proof: 

d 3 

6p 3 = ^(2H 2 A 3 + fc B T m tanh- 1 A) 



= 12H 2 + 2k B T r 

= k B T m 
48k B T, 



4A 2 + (1 - A 2 ) 



2^3 



(1-A 2 ) 



A 



A 3 \ 1 - A 2 
x*dx 



A 3 



(l-x 2 )" 



tanh 1 A 



>0, 



3A 2 



1 



a-A 2 r 



(59) 



where A = (L) T . ■ 

This immediately justitfies the representation of 7r£ ~ 
V2aHl{2p 3 k B T m )^ 1 > 0. Both the diagram in FIG. ||b) 
and analytical expression here show that 



lim 



P3 



(L) Tm ^o k B T„ 



= 0. 



Hence, 



lim 7r£ = 00. 



(60) 



(61) 



This reminds us of the infinite correlation length that 
characterizes the critical point of a continuous phase 
transition - now a limit case of a sequence of first-order 
SLPTs. 

Fortunately, the effectiveness of most of the inferences 
in our model depends on the sign of the parameters in- 
stead of their exact values, so the related physical picture 
of melting will not change even if parameter estimates are 
modified by better methods. 



V. CONCLUSIONS 

In Paper II, we have provided the detailed arguments 
involved in the SLPT model for a 3D fee elemental crys- 
tal. We have used two independent approaches to obtain 
the same equation for the evolution of defect concentra- 
tion. The solution of the model shows that the phe- 
nomenological concerns in Lindemann and Born criteria: 
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atom displacement and lattice softening alone are only 
conducive to, but not decisive in the melting process. 
The two criteria bespeak the underlying cooperative cre- 
ation of interstitials and inter-defect correlations, which 
are crucial for the feasibility of a first-order SLPT marked 
by catastrophes of both rigidity and displacement. The 
sinusoidal correlation waves give rise to interstitial clus- 
ters near the melting point and trigger the heterogeneous 
formation of the SDIs - the pre-liquid droplets satisfying 
Born and Lindemann criteria simultaneously. 

APPENDIX A: THE DEFECT 
CONCENTRATION EQUATION REVISITED 

1. "Chemical Equilibrium" Without H2 and T 



few paragraphs. 

We set out to evaluate the configurational energy dif- 
ference between defective and non-defective cells when 
(1 + |(L) T |)/2 = P and (1 - |(L) T |)/2 = q. Here, P 
is the probability that a site is occupied "correctly" and 
q = 1 — P. We use the following diagrams to aid calcula- 
tion, bearing in mind that each edge in the cube below is 
shared by four cubes and each diagonal line on the sur- 
face of the cube is shared by two cubes. Every vertex of 
the cube (shared by eight cubes) is labeled by the site 
occupancy - the probability that the corresponding site 
is occupied by an atom. 

If NNN correlation is negligible, we may use the fol- 
lowing graph counting to evaluate the average energy 
difference between a non-defective and a defective cell: 



In this subsection, we will provide an alterna- ^- ^ 

tive derivation of the "chemical equilibrium" condition — — = - (£ n0 ii-defective cell - ^defective cell) (Al) 
Eqn. H38|) so as to double-check Eqn. (|23|) . For simplic- 
ity, vibration is provisionally neglected in the following 




(-J 1 +2J 2 )(P-q) 
(-Jx + 2J 2 )|(L)J. 



(A2) 



This value is equal to — (l/4)Ae because for each cube, 
the state of 2 x (1/8) = 1/4 atoms are definite, marked by 
occupancy 1 and respectively. Therefore, the dynamic 



equilibrium between the two types of cell requires that 



P 1 + \{L) T \ 



exp 



(A3) 



where Hi = (3/2) (J a - 2J 2 ). This agrees with Eqs. (H 
and l|38l) when terms associated with H2, T are neglected. 
As \(L) T \ -> 1, Eqn. ijA"3|l implies q ~ exp (-u/2k B T), 
which agrees with the well-established theory of Frenkel 
defects. Here, u = 4Hi is the excitation energy of an 
interstitial defect. 
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2. "Chemical Equilibrium" involving TL2 and T 

When \(L) T \ is sufficiently far from unity, the coop- 
eration and correlation between NNN atoms could not 
be neglected. We can estimate such NNN correlation by 



postulating that the three sites that are NNN neighbor 
to the probability-one-occupied site tend to be occupied 
or be unoccupied simultaneously. Considering this coop- 
eration effect, we should modify —(1/4) As by adding the 
contribution from the following diagrams: 



V 1 



\1 



( 

p 



i ./ 



'3J 2 



P 6 _ 3 ^ p < 
2 2 



3J, 



(P-q)Pq 



^\(L) T \(\(L) T \ 2 



1). 



(A4) 



Accordingly, the chemical equilibrium condition now 
reads: 



1 ~ K L )tI = e - 1 ^ T [2H 1 \{L) T \-iH2\(L) T \{\{L) T \ 2 -l)] 

1 + IWrl 



(A5) 



where Hi = (3/8) J2, in perfect agreement with Eqs. 
(I23[) and l|38|) leaving alone the vibration term. 

The contribution from vibrations could be easily 
considered by the transformation 47ii l_ > — 
ik B T\og{{vi) I {vi)) because averagely speaking, exciting 
one Frenkel pair changes one eigenfrequency from (yi) to 

The chemical equilibrium between defective and non- 
defective cells is reached by stochastic motion of inter- 
stitials, which guarantees sufficient randomness of NN 
pair states (a corollary of short relaxation time of NN 
atom pairs) thereby supporting the former argument that 
holds NN correlation to be negligible. In the mean time, 
NNN atom pairs have considerably long relaxation time 
and their correlation will lead to spatial inhomogeneity, 
which is to be elaborated in the Appendix B. 



APPENDIX B: THE SPATIAL CORRELATION 
OF EXCITATIONS IN DETAIL 

1. The Green Function 

In this subsection, we establish the equation of order 
correlation functions and obtain expressions for the cor- 
relation length at low temperatures. We now start to 
work with the unit system where a — y/2. In this spe- 
cial system, J d 3 r could be interpreted as the summation 
over the 2N sites as well as the volume integration. 



The equation of L (r) spatial correlation Green func- 
tion is obtained by the following procedure [27| : 

(1) We introduce a phenomenological external field 
H (r) and rewrite the free energy functional as 

F [L (r) , H (r)] = F[L (r)] - J d 3 r H (r) L (r) . (Bl) 

(2) We use functional derivative 8 to verify 

S 2 F[L(r),H(r)} 



k B T 
1 

1 



SH (r) 5H (r') 
1 6 2 Q 



1 SQ 1 



Q SH (r) SH (r') Q SH (r) Q SH (r') 
[{L{v)L {v'))-{L{v)) (L(r'))\ 

G(r,r') 



(B2) 



where Q = exp (—F/ksT). 

(3) On the other hand, SF [L (r) , H (r)] /SL (r) = 
implies that 



- 7 V 2 i (r) — p(L (r)) -H(r)= 0. 
Recalling that 

S (L (r)) 5 2 F [L (r) , H (r)] 



SH (t>) 



SH (r) SH (r') 



and 



SH (r') 

we apply the functional differentiation to obtain: 
° = ^)hV 2 iW-p(i(r))^W] 



(B3) 



(B4) 



(B5) 



1 



k R T 



-7V 2 



dp(L(r)) 



dL (r) 



G (r, r') - 5 (r - r') • 
(B6) 
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K, 



K, 



K, 



(a) 




K i >r\ K 



is then 




K 4 K 3 



(b) 



(c) 



FIG. 5: (a) K1+K2+K3 = where |Ki| = |K 2 | = |K 3 |. (b) 
Ki + K 2 +K 3 +K 4 = where |Ki| = |K 2 | = |K 3 | = |K 4 |. 
In general, such "off-plane" term p (Ki) p (K2) p (K3) p (K4) 
does not contain |/i(Ki)| 2 explicitly, (c) "In-plane" four 
wave vectors. In such case, p (Ki) p (Ka) p (K3) p (K4) = 

MRi)! 2 !^)! 2 . 




' » 


ImA 

O 1 




ReA 



(b) 



FIG. 6: (a) This plots the poles (x) and the integration con- 
tour (thick lines) for M e (K) in the complex K plane. Accord- 
ing to Jordan's lemma, only the residues at the poles in the 
upper plane are taken into account, (b) This plots the poles 
(x) and the integration contour (thick lines) for M{K). Note 
that the 4 poles of M £ (K) near the real axis are merged as 
4/2=2 poles in M (K), so every pole on the real axis is wound 
around only half a circle in the contour designed for M (K). 



which is equivalent to Eqn. I|4U|). It is clear that G(r, r') 
is the propagator of Eqn. I|23|). 

For T < T m , where T m is defined by Eqn. 
p(L(r)) is almost linear in the vicinity of (L) T , so 
—dp(L(r)) /dL(r) is almost a constant pi even when 
fluctuation is considered. The solution of the propagator 



G(r,r') 



k B T 



e r/ ^,r = |r - r'| 



(B7) 



— 0, the linearity of p(L(r)) 
T is lost and the corresponding 



However, for T — > T, 
in the vicinity of (L) 
non-Gaussian fluctuations of L (r) in different magni 
tudes have different propagation modes and correlation 
length. A fluctuation characterized by |I/(r)| > \ {L) T | 
decays exponentially (because dp(L(r)) /dL(r) < in 
this case), but fluctuations with |£ (r)| < \(L)\ T could 
be propagated in the form of a sinusoidal wave (because 
dp (L (r)) /dL (r) > in that case). The mean wave- 
length will be estimated in the next subsection. 

Such fluctuation propagation is a consequence of 
plasma instability at T m , and gives rise to collaborative 
atom clusters, the shape of which is nearly spherical, for 
the reason outlined in the next subsection. 



2. The Correlation Functions When T -> T m - 

The mean wavelength of L (r) fluctuation correlation 
could be estimated with the following steps: 

Step 1: Define the Fourier expansion of p (r) as 



(B8) 



K 



so that the complex conjugate of p (K) is p (— K) and 
that 



V M (r) = ^K M (K)e* K - r . 



(B9) 



K 



Expand the free energy functional in the power series of 
p (K) in a volume V = 2N which contains 2N sites: 



F-F [{L) T J 
= |d 3 r { |m 3 (r) p 2 + i^ 4 (r) p 3 + | [V M (r)] 2 } 

= / d3r 5E( !K )( !K 'W K )^ K ') e ' (K+K)r + f E M(K 1 )A i (K 2 ) A1 (K3)e^+ K ^ K 3)- 

[ K,K' Ki,K 2 ,K 3 

+ ? E * ( K i) ( K ^) M (K3) P (K 4 ) el (K 1+ K 2+ K 3+ K 4 ). r 1 

Ki,K 2 ,K 3 ,K4 J 

= y^^ 2 | M (K)| 2 + ^ p(K 1 )p(K 2 )p(K 3 ) + ^ J2 M (Ki) p (K 2 ) p (K 3 ) p (K 4 ) 

K K!+K 2 +K 3 =0 K!+K2+K3+K 4 =0 

(BIO) 
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where K = IKI 



Step 2: Evaluate (|m(K)| ) by 

V(K)| 2 ) = j\ii{K)\ 2 e-^J[^{K) J, 



(Bll) 



K 



Based on the assumption that the fluctuation wave has a definite wave length 2tt£, we can assert that only |K| ~ l/£ 
terms dominate the sum in Eqn. (|B10l) (FIG. Eland Ref. [H) and carry out the summation only on equilateral polygons 
of K. This results in 



F-F[{L) T J=Vj2iK 2 \v(K)\ 2 + 



E Im(k)i 2 i^(k') 



(B12) 



K 



lK|=|K'i«e-! 



where "off-plane" contributions of /x(Ki) p(K.2) /J, (K3) fj, (K4) (FIG. 0(b)) are assumed to be cancelled by phase 
randomness. By using the density of states in K-space, one may find 



E im k ') 



V" 



(27r) ie-'-AK 



|/i(K')r47rX 2 dK 



3tt 2 



when the "spectral width" reads AK — c£ 1 

(ii) 



|K|=|K'K 

which leads to the approximation: 



E iM(K)i 2 iM(Koi 2 ^ yc 3 f t 3) E i^(K)i 4 , 



(B13) 



(B14) 



J d/x (K) 



exp ■ 



d/i (K) exp 



k B T„ 

V 

k B T m 



V|MK)| 2 + 0(M 3 )^ 3yC(c2 + 3) ' — 



4 3vr 2 £ 3 



l/4K)| 4 



t^ 2 Im(k)| 2 



^ 4^ 3^ IMK)I +Q(M) 



( 2k B T m \ 1/2 _ V P3 Vc(c 2 + 3) 3 ( 2k B T m \ 
\ VjK 2 J Ak B T m 3tt 2 £ 3 4 ^ V-yK 2 J 



5/2 



1 

if7 



(B15) 



y d/i(K) |/i(K)| 2 exp- 



r 



k B T„ 



^ 2 HK)| 2 + (M 3 ) + ? ^g_^) 



Im(k)| 4 



d M (K) |^(K)| 2 + 0(^ 



exp 



r 



k B T n 



^ 2 iM(K)r 



1 / 2k B T m 



2 V 



3/2 



l + O 



X 4 



(B16) 



for short range correlations where the K 2 /z 2 term overwhelms the /x 4 term, and this finally leads to 



/I (KM 2 ) = 1 / 2fc B T m \ 3/2 r/ 2fc B T m \ 1/2 Fp 3 ^c(c 2 +3) 3 / 2fc B r m \ 5/2 f I 

jl ; 2VF 7 #V ^v^x 2 ; 4fc B r m 3tt 2 £ 3 4 ^y 7 x 2 y \k 7 



1 y 7 K 2 

2 2k B T rn 

k B T ml K 2 
V 



f V 1 K 2 \ 



F 2 p 3 c(c 2 + 3) 



7 2 if 4 



16fc B T m ^ 2 £ 3 
1 p 3 k B T m _ , 2 



if 2 



4 7T 2 ^ 3 



c ( c2 + 3 )+°(^ 



(B17) 
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Let 



M £ (K) = jK 2 
M (K) = -fK 2 



2^4 1 Pzk B T m , 2 s „ / 1 



7 ^ 4 _lP3M, c(c2 + 3) 



4 7T 2 ^ 3 



(B18) 
(B19) 



where K is a complex variable in general. M e (K) has six poles 33]. Three of these poles lie above the real axis, and 
the other three conjugate poles lie below the real axis. M (K), which is an approximation of M e (K), has four poles 
(FIG. ©). Among the four, there are two poles on the real K axis corresponding to the four poles of M e (K) that lie 
near the real axis. Each of the two poles should be regarded as half a pole that corresponds to the contribution from 
one upper-plane pole in M e (K) when applying the residue theorem to evaluate the contour integral. Therefore, the 
coefficient before the summation of the residues at these two poles is ni instead of 2ni in the equation below. The two 
real-valued poles of M (K) contribute to a sinusoidal wave. The Green function corresponding to this non-Gaussian 
fluctuation is then 



G(r) 



V 



(27T)' 

V 



J d 3 K(MK)|V K ' : 



+ 00 



K 2 dK I , 
'o 



iKr cos ( 



(270* Jo 

I rp p-\-OC 

/ M ;A : A .-in KrdK 
Jo 



2-rr 

sin0d6> / d<p (|/*(K)| S 
o 



2n 2 r 
k B T rn d 



47r 2 r 


dr 




d 


47r 2 r 


Or 




d 



47r 2 r dr 
k R T n 



/+oo 
M e (K) cos KrdK 
-OO 

2iri res(M £ (K) e lKr , K) 

ImK>0 

2ni ^ ies(M (K)e lKr , K) + ni ^ rcs(M (K) e lKr , K) 



= fc-B-^m l e -K r 

8TTjr 



ImK>0 

cos Kor) 



Im/C=0 



(B20) 



Here. 



Therefore, 



K = -t 



which means 



1 p 3 k B T„, 
4 7r 2 7 2 £ 3 

4tt 2 7 2 2 3 / 2 



c(c 2 + 3) 



2V27T 2 7 2 

pzk B T m c(c 2 + 3) a 3 _ pzk B T m o? 



> 



(B21) 



where the NNN atom distance a is restored. (The in- 
equality results from the fact that c < 1 - that is, the 
maximum spectral width cannot exceed Kq.) 



2V27T 3 7 2 
P3k B T m a 3 



(B22) 



gives an estimate of the minimal diameter of unstable 
clusters in the solid when T = T m . 
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